用LeafCutter探索转录组数据的可变剪切

该软件早在2016年就公布了，发表在biorxiv预印本上面，但直到2017年的双11，才发表在NG上面，文章是 : Annotation-free quantification of RNA splicing using LeafCutter 最大的特点应该是不需要参考基因组的基因注释信息了吧，就是gtf/gff文件可以省略，当然，比对还是需要的。它还有另外一个非常重要的功能，splicing quantitative trait loci (sQTLs) 但是跟我目前关系不大， 就不介绍了。

本教程，首发于 生信菜鸟团博客：http://www.bio-info-trainee.com/2949.html

背景介绍

目前主流的探究转录组数据的可变剪切的算法要么是基于estimate isoform ratios 或者 exon inclusion levels ，但是挑战还是蛮多的，可变剪切本跟正常转录本重合的比例很大，技术误差也是有的，依赖于基因现有的注释信息，既不准确，也不完全。所以作者开发了LeafCutter工具。

LeafCutter workflow.

• First, short reads are mapped to the genome. When SNP data are available, WASP should be used to filter allele-specific reads that map with a bias.

• Next, LeafCutter extracts junction reads from.bam files, identifies alternatively excised intron clusters, and summarizes intron usage as counts or proportions.

• Finally, LeafCutter identifies intron clusters with differentially excised introns between two user-defined groups by using a Dirichlet-multinomial model, or maps genetic variants associated with intron excision levels by using a linear model.

作者在Genotype-Tissue Expression (GTEx) Consortium数据集上面测试了，并且把结果跟 GENCODE v19, Ensembl, and UCSC 着3大主流的基因注释信息数据库比较。还在其它数据库里面验证了，数据下载地址是：dbGaP under accession phs000424.v6.p1 (GTEx), GEO under accession GSE41637 (RNA-seq data from mammalian organs), and ENA under accession PRJEB3366 (Geuvadis).

软件下载地址：

• LeafCutter software, https://github.com/davidaknowles/leafcutter;

• LeafViz visualizations, https://leafcutter.shinyapps.io/leafviz/;

• rheumatoid arthritis summary statistics, http://plaza.umin.ac.jp/yokada/datasource/software.htm.

软件安装及使用

最简单的就是conda进行安装了：

1. conda install -c davidaknowles r-leafcutter

如果安装失败，可能需要单独为它创建一个环境。

不过，它本身就是一个R包，所以在个人电脑里面的rstudio里面安装即可。

1. if (!require("devtools")) install.packages("devtools", repos='http://cran.us.r-project.org')

2. devtools::install_github("davidaknowles/leafcutter/leafcutter")

但是源代码里面有一些脚本和测试数据，所以还是要下载看看

1. mkdir -p ~/biosoft

2. cd ~/biosoft

3. git clone https://github.com/davidaknowles/leafcutter

4. cd leafcutter

5. ## 需要修改里面的一个脚本 scripts/bam2junc.sh 把软件路径增添进去即可

里面又是perl又是python的，感觉他们团队开发环境不统一。

第一步:bam2junc

比对一般来说，优先选择STAR等支持跨越内含子的转录组比对工具得到bam文件，运行下面的脚本即可进行批量转换：

1. cat bam_path.txt |while read id

2. do

3. file=$(basename $id )

4. sample=${file%%.*}

5.    echo Converting $id to $sample.junc

6.    sh /public/biosoft/leafcutter/scripts/bam2junc.sh  $id $sample.junc

7. done

得到的junc文件如下:

1. chr7    134840725   134843893   .   1   -

2. chr2    234355442   234355737   .   1   +

3. chr4    37828435    37831585    .   13  +

4. chr19    39101772    39101882    .   5   +

5. chr11    109735445   109827551   .   19  +

6. chr18    48458730    48465939    .   8   -

7. chr12    82751048    82752457    .   12  -

8. chr15    51018323    51018517    .   14  -

9. chr1    247323115   247335149   .   2   +

10. chr10    92920631    92982445    .   1   +

这个步骤有点耗时，所有的junc文件地址需要保存给下一步使用

第二步：Intron clustering

这个步骤，需要python2.7版本，这个是python的一个大坑，到现在版本仍然不统一。

1. ls *.junc >test_juncfiles.txt

2. python /public/biosoft/leafcutter/clustering/leafcutter_cluster.py -j test_juncfiles.txt -m 50 -o testYRIvsEU -l 500000

几分钟就运行完毕。

得到的比较重要的文件如下：

1. 1.3M Jan  4 17:45 testYRIvsEU_perind.counts.gz

2. 680K Jan  4 17:45 testYRIvsEU_perind_numers.counts.gz

3. 5.0M Jan  4 17:45 testYRIvsEU_pooled

4. 540K Jan  4 17:45 testYRIvsEU_refined

5. 877 Jan  4 17:45 testYRIvsEU_sortedlibs

6. 854 Jan  4 17:43 test_juncfiles.txt

值得注意的是 testYRIvsEU_perind_numers.counts.gz 文件，里面每一行都是一个内含子，每一列都是一个样本，写明了它们的表达值，这些数值就可以用来做可变剪切分析。

1. #  zcat testYRIvsEU_perind_numers.counts.gz |tail

2. chr8:145651155:145651305:clu_6538 21 14 19 8 9 0 13 33 0 0 4 0 5 8 12 0 12 34 15 0 0 10 11

3. chr8:145651155:145651409:clu_6538 1021 611 186 190 294 284 681 89 222 57 257 363 694 807 523 44 469 812 926 71 80 260 214

4. chr8:145652362:145653872:clu_6539 1265 694 132 74 302 71 178 34 44 12 63 122 230 218 472 6 146 1421 1084 16 14 83 46

5. chr8:145652654:145653872:clu_6539 48 24 56 0 26 0 13 0 2 5 2 0 3 19 17 0 2 8 64 0 0 3 0

6. chr8:145652674:145653872:clu_6539 18 26 0 0 0 7 2 0 5 0 0 0 1 6 11 0 3 34 37 0 0 9 6

7. chr8:146017525:146017630:clu_6540 2 3 44 0 2 12 4 0 0 0 22 5 9 10 2 0 1 9 11 0 0 1 0

8. chr8:146017525:146017751:clu_6540 1067 671 620 41 295 347 224 89 62 33 262 136 229 223 356 17 288 480 1842 9 35 70 23

9. chr8:146076780:146078224:clu_6541 18 3 0 0 17 17 8 0 0 3 2 3 16 6 12 0 4 45 29 9 0 10 2

10. chr8:146076780:146078378:clu_6541 22 17 0 0 0 3 1 0 0 0 3 2 15 7 2 0 7 62 55 0 0 4 0

11. chr8:146076780:146078757:clu_6541 10 1 16 0 12 52 0 0 11 0 24 9 27 3 0 0 7 0 28 0 0 2 0

第三步：制作分组矩阵进行差异分析

避免暴露我真实的项目，这里就给作者的示例文件吧：

1. RNA.NA18486_YRI.chr1.bam YRI

2. RNA.NA18487_YRI.chr1.bam YRI

3. RNA.NA18488_YRI.chr1.bam YRI

4. RNA.NA18489_YRI.chr1.bam YRI

5. RNA.NA18498_YRI.chr1.bam YRI

6. RNA.NA06984_CEU.chr1.bam CEU

7. RNA.NA06985_CEU.chr1.bam CEU

8. RNA.NA06986_CEU.chr1.bam CEU

9. RNA.NA06989_CEU.chr1.bam CEU

10. RNA.NA06994_CEU.chr1.bam CEU

很简单的两列文件，说明每一个样本属于哪个组即可。

1. /public/biosoft/leafcutter/scripts/leafcutter_ds.R --num_threads 4 \

2. --exon_file=/public/biosoft/leafcutter/leafcutter/data/gencode19_exons.txt.gz \

3. testYRIvsEU_perind_numers.counts.gz group_info.txt

这里的 group_info.txt 就是自己制作好的分组矩阵。值得提醒的是，上面的文件有且只能有2个分组，这样软件才知道怎么样去比较，如果自己的分组很多，可以考虑制作多个分组文件，运行多次。

当然，上面的脚本已经没有必要在linux服务器里面运行啦。

既然有了内含子的表达矩阵，又有了分组信息，差异分析根本就不会消耗多少计算资源，全部下载到自己的电脑里面去做吧。

自己打开文件 /public/biosoft/leafcutter/scripts/leafcutter_ds.R 就明白了整个流程。

也是几分钟就完成了全部结果。

1. Running differential splicing analysis...

2. Differential splicing summary:

3.                                             statuses Freq

4. 1 <2 introns used in >=min_samples_per_intron samples  425

5. 2                          <=1 sample with coverage>0   62

6. 3               <=1 sample with coverage>min_coverage  939

7. 4                            Not enough valid samples 3047

8. 5                                             Success 2068

9. Saving results...

10. Loading exons from /Users/jmzeng/biosoft/leafcutter/leafcutter/data/gencode19_exons.txt.gz

11. All done, exiting

得到的文件里面，需要详细了解的是 leafcutter_ds_cluster_significance.txt 主要靠自己看readme啦。

第四步：可视化那些可变剪切

也是包装好的脚本。

1. /Users/jmzeng/biosoft/leafcutter/scripts/ds_plots.R -e  /Users/jmzeng/biosoft/leafcutter/leafcutter/data/gencode19_exons.txt.gz testYRIvsEU_perind_numers.counts.gz   group_info.txt leafcutter_ds_cluster_significance.txt -f 0.05

所有的可变剪切形式都会可视化在一张图里面。

可视化函数还是有很多参数是可以自己调整的，值得学习

leafcutter得到剪切形式

上面的剪切形式略微有点复杂。

leafcutter得到剪切形式

这个稍微简单一些，可以看到虽然剪切形式没有变化，但是不同的转录本表达量差异很大。

leafcutter得到剪切形式

最后这个有可变剪切形式。